Astronomical pacing of the Lau Biogeochemical Event captured in the Kosov Quarry section
A study on the imprint of astronomical cycles in the Kosov quarry spanning the Lau Event during the Silurian
Authors
Michiel Arts1
Anne-Christine da Silva1
1 - SEDICLIM lab, Department of Geology, University of Liege, Belgium
Published
October 1, 2026
Note
This study is based onthe induration record extracted from the litholog of Fryda, Jiri & Lehnert, Oliver & Frýdová, Barbora & Farkas, Juraj & Kubajko, Michal. (2020). Carbon and sulfur cycling during the mid-Ludfordian anomaly and the linkage with the late Silurian Lau/Kozlowskii Bioevent. Palaeogeography, Palaeoclimatology, Palaeoecology. 564. doi:10.1016/j.palaeo.2020.110152
The Kosov Quarry section in the Prague Basin preserves one of the most complete Silurian successions spanning the Lau Biogeochemical Event. This interval records major perturbations in the carbon cycle, climate cooling, redox fluctuations, and faunal turnover, yet the pacing and duration of these changes remain poorly constrained. Leveraging its exceptional completeness, we developed a new astrochronology based on the induration record, which allows the identification of cycles ranging from the half-precession to the 405-kyr eccentricity cycle. The anchored astrochronology that indicates that the Lau Biogeochemical Event started at 424.06 ± 0.55 Ma and lasted 0.87 ± 0.09 Myr, while the associated Lau-Kozlowski Bio-event was much shorter at 30 ± 10 kyr. Astronomical forcing is expressed not only in the induration record but also in the rate of change (‰/kyr) of the δ13Ccarb curve and in the lag-1 sea-level proxy. The δ13Ccarb rate-of-change record shows a strong 100-kyr eccentricity imprint during the onset of the Lau Event, reaching a maximum of 0.09 ‰/kyr. Similarly, the lag-1 sea-level record is dominated by the 405-kyr eccentricity cycle and appears in-phase with insolation, consistent with (405-kyr) eccentricity-paced glacio-eustatic sea-level fluctuations during the Ludlow.
1. Introduction
The Lau Biogeochemical Event is associated with one of the largest carbon cycle perturbations of the entire Phanerozoic. The Lau Biogeochemical Event is marked by a large positive δ¹³C excursion, faunal turnover, and evidence for sea-level and redox fluctuations (Frýda and Manda, 2013; Frýda et al., 2020). These changes indicate a tight coupling between the carbon cycle, climate, and oceanographic processes; however, the tempo and (astronomical) pacing of the Event remain poorly constrained. Without a precise temporal framework, it is not possible to determine whether the Lau reflects a brief perturbation or a protracted interval of instability.
To address this, we will conduct a cyclostratigraphic study on the Kosov quarry section, which provides a nearly complete archive, allowing us to establish a high-resolution astrochronological framework for the Lau Biogeochemical Event. The induration record of Frýda et al. (2020) will form the basis for constructing the astrochronology, which will enable us to assign durations and associated uncertainties to the Lau Biogeochemical Event and its associated (sub)intervals. The lag-1 sea-level proxy will be derived from the tuned induration record, which will provide constraints and insights into sea-level variability during the Lau Biogeochemical event. The rate of change (‰/kyr) calculated from the δ13Ccarb record will be used to gain insights into carbon cycle dynamics during the Lau Biogeochemical Event. Furthermore, the imprint of astronomical cycles on the δ13Ccarb rate of change (‰/kyr) and lag-1 sea-level will be delineated to understand the role that astronomical cycles might have played in pacing climate, organic carbon burial, and sea-level variations during the Lau Biogeochemical Event.
1.1 Geological setting
Our study focuses on upper Silurian (Ludfordian) strata of the Prague Basin, a sedimentary succession within the Barrandian (Teplá-Barrandian terrane) (Frýda and Manda, 2013; Frýda et al., 2020, 2021a) (Figure 1). During the Ludfordian, this block was located at ~25–30°S (Tasáryová et al., 2014; Scotese, 2021). Facies distribution indicates hemipelagic deposition surrounded by shallow- to deeper-marine environments (HORNÝ, 1955; Kříž, 1991). It remains debated whether the Barrandian’s geodynamic setting belongs to the peri-Gondwanan terrains or is its own isolated Perunica microcontinent (Fatka and Mergl, 2009).
Code
par(mar =c(4, 2, 2, 2))# Empty plot to define margins and space in the panelplot(1, type ="n", axes =FALSE, xlab ="", ylab ="", main="A")# Get the coordinates of the current plot regionusr <-par("usr") # c(x1, x2, y1, y2)# Draw the image to fill the entire panelrasterImage(img_2, usr[1], usr[3], usr[2], usr[4])
Figure 1. Regional overview of the Kosov section. Copied from Fryda et al. (2020)
Code
par(mar =c(4, 6, 2, 2))
The studied Kosov section (No.JF195; 49°56′12.2″N, 14°03′20.1″E) provides one of the most complete global records of the mid-Ludfordian carbon isotope excursion (Frýda and Manda, 2013; Frýda et al., 2020, 2021a) (Figure 2).The succession comprises carbonates with marly interbeds and shales, including micritic mudstones, skeletal limestones, and crinoidal grainstones (Kříž, 1992; Lehnert et al., 2007). The site has been widely studied for its sedimentology, palaeontology, and faunal changes across the Lau/Kozlowskii Bioevent (LKB) and Mid-Ludfordian Carbon Isotope Excursion (MLCIE) (Kříž, 1992; Štorch, 1995; Manda and Kříž, 2006; Frýda and Manda, 2013).
Figure 2. The original litholog, chronostratigraphy, event zonation, induration record and the δ13Ccarb curve. A. The Litholog of Frýda et al. (2020). B. Chronostratigraphic subdivision. C. Event stratigraphic intervals. D. The induration record extracted from the litholog E. The δ13Ccarb record of Frýda et al. (2020)
2. Materials and Methods
A cyclostratigraphic analysis was conducted on the induration record extracted from the litholog published in Frýda et al. (2020) (see https://github.com/stratigraphy/Kosov-cyclostrat for the R code). The induration record is considered as a proxy for lithification intensity and sedimentary competence, with low values corresponding to soft shales and marl and high values corresponding to competent limestone beds. The analyses leveraged the functions of the WaverideR (Arts, 2023; Arts et al., 2024) and astrochron R packages (Meyers, 2015, 2019) to perform the cyclostratigraphic study.
The cyclostratigraphic analysis began by applying the Time Optimisation (TimeOpt) function of the astrochron R package to the induration record, yielding a first-order estimate of the sedimentation rate (Meyers, 2015). The TimeOpt function was run twice: once to evaluate the 405-kyr eccentricity amplitude modulation of the 100-kyr eccentricity band and the concentration of power at the 100-kyr and 405-kyr eccentricity frequencies, and once to evaluate the eccentricity amplitude modulation of the precession band and the concentration of power at the precession and eccentricity frequencies. The optimal fit of the astronomical cycles was evaluated for a range of sediment accumulation rates spanning between 0.1 and 4 cm/kyr, which covers the minimum and maximum estimates for the durations of the studied interval (Cramer et al., 2015; Melchin et al., 2020; Arts et al., 2024).
The first-order sedimentation rate estimate from the TimeOpt runs was used to extract the 405-kyr eccentricity cycle from the induration record. The extracted cycle was then used to conduct a minimal tuning, where the distance between two successive peaks of the extracted cycle was set to the duration of the interpreted astronomical cycle. This duration was subsequently converted into an age model and a sedimentation-rate curve. The sedimentation-rate curve was then used to convert cm/kyr to the period (m) of different astronomical cycles. These cycles were then plotted on top of the CWT scalogram to verify whether the curves passed through areas of high spectral power, validating the presence of astronomical cycles. When this was the case, the period (m) of the 405-kyr eccentricity cycle was tracked in the CWT.
The tracked period (m) curve of the 405-kyr eccentricity cycle was then used to create a floating astrochronology with an assigned uncertainty based on the analytical uncertainty of the Wavelet (Arts et al., 2024). This floating age model with uncertainty was then anchored to the astrochronologically calibrated age of the Ludlow–Pridoli boundary (423.03 ±0.53 Ma) (Arts et al., 2024). The anchored numerical age model was then used to assign ages and durations (with uncertainty) to the record of the Kosov quarry section and its subsections.
The mean age model was then used to convert the induration, allowing for the study of the imprint of astronomical cycles in the induration record in the time domain. Special attention was paid to the phase relationship between the 100-kyr eccentricity cycle directly extracted from the record and the 100-kyr eccentricity cycle extracted from the Hilbert transform of the precession cycle, and between the 405-kyr eccentricity cycle directly extracted from the record and the 405-kyr eccentricity cycle extracted from the Hilbert transform of the 100-kyr eccentricity cycle. The phase relationship between these cycles allowed the inference of the phase relationship between cycles extracted from the proxy record and the true phase of the astronomical cycle, as well as its relationship with changes in insolation (Hinnov, 2000; Laskar et al., 2004).
The δ13Ccarb record was also converted to the time domain which allowed the for the calculation of the rate of change (‰/kyr) in the tuned δ13Ccarb record. This record was then further investigated for the imprint of astronomical cycles.
A lag-1 curve was also generated based by applying a windowed lag-1 autocorrelation coefficient Monte Carlo simulation on the induration record in the time domain. This record functions as a proxy record for the sea-level curve (Li et al., 2018). The imprint of astronomical cycles in the lag-1 record was investigated and compared to the imprint of astronomical cycles in the induration record to understand the relationship between astronomical forcing and sea-level change during the Lau Biogeochemical Event.
3. Results
3.1.TimeOpt-derived sedimentation rates
The TimeOpt function was run twice, once to optimize the 405-kyr eccentricity amplitude modulation of the 100-kyr eccentricity band, and the concentration of power at the 100-kyr and 405-kyr eccentricity frequencies and once to optimize the eccentricity amplitude modulation of the precession band, and the concentration of power at the precession and eccentricity frequencies. The TimeOpt results indicate that the optimal sedimentation rate is around 3 cm/ka, with the peak spanning sedimentation rates between 2 and 4 cm/ka (see Figure 3).
Code
#run timeOpt#we need to run it twice for a correct plot later ontargetP_AM=c(405.7, 130.7, 123.8, 98.9,94.9)targetP=c(21,17)res1=timeOpt(dat = induration,sedmin =0.1,sedmax =4,numsed =500,fit =2,fitModPwr =TRUE,targetE=targetP_AM,roll =10^4,output=1)res2 <-timeOpt(dat = induration,sedmin =0.1,sedmax =4,numsed =500,fit =2,fitModPwr =TRUE,targetE=targetP_AM,roll =10^4,output=2,genplot=TRUE)res3=timeOpt(dat = induration,sedmin =0.1,sedmax =4,numsed =500,fit =1,fitModPwr =TRUE,targetP=targetP,roll =10^4,output=1)res4 <-timeOpt(dat = induration,sedmin =0.1,sedmax =4,numsed =500,fit =1,fitModPwr =TRUE,targetP=targetP,roll =10^4,output=2,genplot=TRUE)
Code
layout(matrix(c(1, 2, 3, 4,5,6,7,8), nrow =4, ncol =2, byrow =TRUE)) par(mar =c(4, 4, 3, 3))# --- Panel A ---plot(res1[, 1], res1[, 2], cex =0.75, cex.lab =1.2, cex.main =1.3, col ="red", xlab ="", ylab ="", main =expression(paste(bold("Fit: "), {"r"^2 }["envelope"], " (red) and ", {"r"^2 }["power"], " (gray)")))axis(2, col.axis ="red")mtext(expression("r"^2), side =2, line =2, cex =0.9)mtext("Sedimentation rate (cm/ka)", side =1, line =2.3, cex =0.8)par(new =TRUE)plot(res1[, 1], res1[, 3], col ="#00000064", xlab ="", ylab ="", type ="l", axes =FALSE, lwd =2)axis(4, ylim =c(0, max(res1[, 3])), lwd =1, col ="black")usr <-par("usr")text(usr[1], usr[4] -0.05*(usr[4]-usr[3]), "A", pos =4, font =2, cex =1.5)# --- Panel B ---plot(res1[, 1], res1[, 4], type ="l", lwd =2, cex.lab =1.2, cex.main =1.3, col ="black", xlab ="", ylab ="", main =expression(paste(bold("Optimal Fit: "), {"r"^2 }["opt"])))mtext(expression({"r"^2}["opt"]), side =2, line =1.9, cex =0.9)mtext("Sedimentation rate (cm/ka)", side =1, line =2.3, cex =0.8)usr <-par("usr")text(usr[1], usr[4] -0.05*(usr[4]-usr[3]), "B", pos =4, font =2, cex =1.5)# --- Panel C ---ylim1 =c(min(res2[,4], res2[,5]), max(res2[,4], res2[,5]))plot(res2[,1], res2[,4], cex =0.5, cex.lab =1.2, cex.main =1.2, xlab ="", ylab ="", main ="Envelope (red); Reconstructed Ecc. Model (black)", col ="red", type ="l", ylim = ylim1)lines(res2[, 1], res2[,5], lwd =1.5)mtext("Std. Value", side =2, line =2, cex =0.9)mtext("Time (ka)", side =1, line =2.3, cex =0.8)usr <-par("usr")text(usr[1], usr[4] -0.05*(usr[4]-usr[3]), "C", pos =4, font =2, cex =1.5)# --- Panel D ---ylim1 =c(min(res2[, 3], -1* res2[, 4]), max(res2[, 3], res2[, 4]))plot(res2[,1],res2[,3], col ="blue", cex =0.5, cex.lab =1.2, cex.main =1.3, xlab ="", ylab ="", main ="Envelope (red); Filtered Data (blue)", ylim = ylim1)lines(res2[,1],res2[,3], col ="blue")lines(res2[,1],res2[,4], col ="red")lines(res2[,1], -1*res2[,4], col ="red")abline(h =0, col ="black", lty =3)mtext("Std. Value", side =2, line =2, cex =0.9)mtext("Time (ka)", side =1, line =2.3, cex =0.8)usr <-par("usr")text(usr[1], usr[4] -0.05*(usr[4]-usr[3]), "D", pos =4, font =2, cex =1.5)# --- Panel E ---plot( res3[, 1], res3[, 2], cex =0.75, cex.lab =1.2, cex.main =1.3, col ="red", xlab ="", ylab ="", main =expression(paste(bold("Fit: "), {"r"^2 }["envelope"], " (red) and ", {"r"^2 }["power"], " (gray)")))axis(2, col.axis ="red")mtext(expression("r"^2), side =2, line =2, cex =0.9)mtext("Sedimentation rate (cm/ka)", side =1, line =2.3, cex =0.8)par(new =TRUE)plot( res3[, 1], res3[, 3], col ="#00000064", xlab ="", ylab ="", type ="l", axes =FALSE, lwd =2)axis(4, ylim =c(0, max( res3[, 3])), lwd =1, col ="black")usr <-par("usr")text(usr[1], usr[4] -0.05*(usr[4]-usr[3]), "E", pos =4, font =2, cex =1.5)# --- Panel F ---plot( res3[, 1], res3[, 4], type ="l", lwd =2, cex.lab =1.2, cex.main =1.3, col ="black", xlab ="", ylab ="", main =expression(paste(bold("Optimal Fit: "), {"r"^2 }["opt"])))mtext(expression({"r"^2}["opt"]), side =2, line =1.9, cex =0.9)mtext("Sedimentation rate (cm/ka)", side =1, line =2.3, cex =0.8)usr <-par("usr")text(usr[1], usr[4] -0.05*(usr[4]-usr[3]), "F", pos =4, font =2, cex =1.5)# --- Panel G ---ylim1 =c(min( res4[,4], res4[,5]), max( res4[,4], res4[,5]))plot( res4[,1], res4[,4], cex =0.5, cex.lab =1.2, cex.main =1.2, xlab ="", ylab ="", main ="Envelope (red); Reconstructed Ecc. Model (black)", col ="red", type ="l", ylim = ylim1)lines( res4[, 1], res4[,5], lwd =1.5)mtext("Std. Value", side =2, line =2, cex =0.9)mtext("Time (ka)", side =1, line =2.3, cex =0.8)usr <-par("usr")text(usr[1], usr[4] -0.05*(usr[4]-usr[3]), "G", pos =4, font =2, cex =1.5)# --- Panel H ---ylim1 =c(min( res4[, 3], -1* res4[, 4]), max( res4[, 3], res4[, 4]))plot( res4[,1], res4[,3], col ="blue", cex =0.5, cex.lab =1.2, cex.main =1.3, xlab ="", ylab ="", main ="Envelope (red); Filtered Data (blue)", ylim = ylim1)lines( res4[,1], res4[,3], col ="blue")lines( res4[,1], res4[,4], col ="red")lines( res4[,1], -1* res4[,4], col ="red")abline(h =0, col ="black", lty =3)mtext("Std. Value", side =2, line =2, cex =0.9)mtext("Time (ka)", side =1, line =2.3, cex =0.8)usr <-par("usr")text(usr[1], usr[4] -0.05*(usr[4]-usr[3]), "H", pos =4, font =2, cex =1.5)
Figure 3. TimeOpt results. A. Fit: r²[envelope] (red) and r²[power] (grey) for TimeOpt run 405 vs 100-kyr ecc. B. Optimal Fit: r²[opt] for TimeOpt run 405 vs 100-kyr ecc. C. Envelope (red); Reconstructed Ecc. Model (black) for TimeOpt run 405 vs 100-kyr ecc. D. Envelope (red); Filtered Data (blue) for TimeOpt run 405 vs 100-kyr ecc. E. Fit: r²[envelope] (red) and r²[power] (grey) for TimeOpt run 100-kyr eccentricity versus precession. F. Optimal Fit: r²[opt] for TimeOpt run 100-kyr eccentricity versus precession. G. Envelope (red); Reconstructed Ecc. Model (black) for TimeOpt run 100-kyr ecc. versus precession. H. Envelope (red); Filtered Data (blue) for TimeOpt run 100-kyr eccentricity versus precession.
3.2 Extract the 405-kyr cycle and generate a sedimentation-rate curve
Based on the TimeOpt results, the 405-kyr eccentricity cycle is extracted from the induration record using a Taner bandpass filter, which spans from 4.86 to 16.2, corresponding to the peak of the optimal sedimentation rate range (1.2-4 cm/ka) as determined by the TimeOpt analysis results. Next, a sedimentation-rate curve is constructed using the minimal tuning technique, in which the distance between two successive peaks of the extracted cycle is set to the duration of the interpreted astronomical cycle. This duration is subsequently used to generate the sedimentation-rate curve (see Figure 4).
Figure 4. Minimal tuning. A. Induration record and the extracted 405-kyr eccentricity cycle. B. Sedimentation rate and sedimentation rate (cm/kyr) based on the minimal tuning.
3.3 Tracking the period (m) of the 405-kyr eccentricity cycle in the the CWT scalogram
The sedimentation rate curve attained from the minimal tuning was then used to convert cm/kyr to period (m) of different astronomical cycles. These cycles were then plotted on top of the CWT scalogram (see figure 4). We can see that the curves plot through areas of high spectral power corresponding to astronomical cycles that should be present according to the TimeOpt-derived sedimentation rate modelling runs. Next, the period (m) of the 405-kyr eccentricity cycle was tracked in the CWT.
Figure 5. The Continuous Wavelet Transform (CWT) of the induration record extracted from the litholog of Frýda et al. (2020) and the period (m) of astronomical cycles based on the sedimentation rates attained from the minimal tuning and the tracked period (m) of the 405-kyr eccentricity cycle. Black dotted = tracked period (m) the 405-kyr eccentricity cycle Red = the 405-kyr eccentricity cycle (minimal tuning based) Purple = the 100-kyr eccentricity cycle (minimal tuning based) Blue = the obliquity cycle (minimal tuning based) Dark green = the precession cycle (minimal tuning based)
3.4 Combine the tracked period (m) and a tie-point curve to create an anchored astrochronoloy
The tracked period (m) curve of the 405-kyr eccentricity cycle was then be used to create a floating astrochronology with an assigned uncertainty based on the analytical uncertainty of the Wavelet (Arts et al., 2024). This floating age model with uncertainty was then anchored to the astrochronologically calibrated age of the Ludlow-Pridoli boundary (423.03 ±0.53 Ma) (Arts et al., 2024). This anchored numerical age model could then be used to assign ages and durations (with uncertainty) to the record of the Kosov quarry section and its subsections
Figure 6. Astrochronological numerival age- depth model.The black line is the mean age-depth curve, and the red and blue lines are time plus and minus two standard deviations (2σ)
3.5 Ages and durations for intervals including uncertainty
The age model provides durations with uncertainty for the chrono, litho-, and event-stratigraphic intervals previously identified in the succession (Frýda and Manda, 2013; Frýda et al., 2021a). The calculated ages, durations, and their associated uncertainties were rounded to the nearest 10 kyr to align with the precision of the anchor point.
Code
kosov_time<-cbind(curve2time_unc_anchor_res[[1]],rowMeans(results))kosov_sd<-rowSds(as.matrix(results))for(i in1:nrow(my_data)){row_nr_1 <- DescTools::Closest(kosov_time[,1],my_data[i,2],which =TRUE)row_nr_2 <- DescTools::Closest(kosov_time[,1],my_data[i,3],which =TRUE)my_data[i,4] <-round(mean(as.matrix(results[row_nr_1,]),na.rm =TRUE)/10,0)/100my_data[i,5] <-round(2*sd(as.matrix(results[row_nr_1,]),na.rm =TRUE)/10,0)/100my_data[i,6]<-round(mean(as.matrix(results[row_nr_2,]),na.rm =TRUE)/10,0)/100my_data[i,7] <-round(2*sd(as.matrix(results[row_nr_2,]),na.rm =TRUE)/10,0)/100my_data[i,8] <-round(mean(as.matrix(na.omit(results[row_nr_1,]-results[row_nr_2,])))/10,0)*10my_data[i,9] <-round(2*sd(results[row_nr_2,]-results[row_nr_1,],na.rm =TRUE)/10,0)*10}for(i in1:nrow(my_data)){for( j in8:ncol(my_data)){if (my_data[i,j]<=10){ my_data[i,j] <-10 }}}bottom_record <- kosov_time[DescTools::Closest(kosov_time[,1],my_data[1,2],which=TRUE),2]base_pridoli <- kosov_time[DescTools::Closest(kosov_time[,1],my_data[8,2],which=TRUE),2]Top_record <- kosov_time[DescTools::Closest(kosov_time[,1],my_data[1,3],which=TRUE),2]base_LAU<- kosov_time[DescTools::Closest(kosov_time[,1],my_data[4,2],which=TRUE),2]base_S_zone<- kosov_time[DescTools::Closest(kosov_time[,1],my_data[5,2],which=TRUE),2]base_F_zone<- kosov_time[DescTools::Closest(kosov_time[,1],my_data[6,2],which=TRUE),2]Top_F_zone<- kosov_time[DescTools::Closest(kosov_time[,1],my_data[6,3],which=TRUE),2]base_event <- kosov_time[DescTools::Closest(kosov_time[,1],my_data[2,2],which=TRUE),2]# Combine uncertainty with age and duration valuesmy_data$age_bottom_combined <-paste0(round(my_data[,4], 3), " ± ", round(my_data[,5], 2))my_data$age_top_combined <-paste0(round(my_data[,6], 3), " ± ", round(my_data[,7], 2))my_data$duration_combined <-paste0(round(my_data[,8], 0), " ± ", round(my_data[,9], 0))# Reorder the data frame BEFORE passing to gtmy_data_out <- my_data |> dplyr::select( Interval, position_bottom, position_top, age_bottom_combined, age_top_combined, duration_combined )# Create gt tabletab <-gt(my_data_out) |>cols_label(Interval ="Interval",position_bottom ="position (m) bottom",position_top ="position (m) top",age_bottom_combined ="Age bottom (Ma) ±2sd ",age_top_combined ="Age top (Ma) ±2sd ",duration_combined ="duration (kyr) ±2sd" ) |>tab_header(title ="Table 1. Ages and durations based on the anchored astrochronological age model" )tab
Table 1. Ages and durations based on the anchored astrochronological age model
Interval
position (m) bottom
position (m) top
Age bottom (Ma) ±2sd
Age top (Ma) ±2sd
duration (kyr) ±2sd
Entire record
-7.2
25.0
424.6 ± 0.56
422.8 ± 0.53
1800 ± 190
Lau Biogeochemical Event
-0.4
20.0
424.06 ± 0.55
423.19 ± 0.53
870 ± 90
LKB Event (trace metal peak)
-0.4
0.0
424.06 ± 0.55
424.03 ± 0.55
30 ± 10
R-zone
0.0
1.4
424.03 ± 0.55
423.94 ± 0.54
90 ± 10
S-Zone
1.4
11.6
423.94 ± 0.54
423.64 ± 0.54
310 ± 30
F-Zone
11.6
20.0
423.64 ± 0.54
423.19 ± 0.53
440 ± 50
Ludlow part section
0.0
22.0
424.03 ± 0.55
423.03 ± 0.53
1000 ± 110
Pridoli part section
22.0
25.0
423.03 ± 0.53
422.8 ± 0.53
230 ± 20
3.6 Extract and plot the record and the astronomical cycles in the time domain
The age model was used to convert the induration and δ13Ccarb record to the (numerical) time domain. Next, the CWT was performed on the tuning induration record. Peaks corresponding to the 8.5-kyr half precession, ~16-kyr precession, 28-kyr obliquity,50-kyr eccentricity, 100-kyr eccentricity, 200-kyr eccentricity and 405-kyr eccentricity cycle can be observed in the CWT scalogram (see Figure 7).Next the precession, obliquity, 100-kyr eccentricity and 405-kyr eccentricity were also extracted from the tuned induration record (see Figure 7). The Hilbert transform was applied to the precession, obliquity and 100-kyr eccentricity cycles extracted from the tuned induration record. The 405-kyr eccentricity cycle was also extracted from the Hilbert transform-derived amplitude record of the 100-kyr eccentricity and compared to the 405-kyr eccentricity directly extracted from the record, showing an anti-phased relationship (see Figure 7).
Figure 7. Induration proxy record in the time domain, including the wavelet scalogram of the induration record and astronomical cycles extracted from said record. A Stages. B. Event zone subdivision C. The δ13Ccarb record. D. The Induration record. E. Wavelet scalogram of the induration record with the average spectral power on top. The black vertical lines in the wavelet scalograms are durations of known astronomical cycles. From left to right, these cycles are the 8.5-kyr half precession, ~16-kyr precession, 28-kyr obliquity,50-kyr eccentricity, 100-kyr eccentricity, 200-kyr eccentricity and the 405-kyr eccentricity cycle. F. The black line is a 405-kyr eccentricity cycle extracted from the Induration record. The red line is the 405-kyr eccentricity cycle extracted from the Hilbert transform of the 100-kyr eccentricity cycle of the Induration record. G. The 100-kyr eccentricity cycle was extracted from the Induration record (black line) and the Hilbert transform of the 100-kyr eccentricity cycle (red line). H.) Obliquity cycle extracted from the Induration record (black line), and the Hilbert transform of the obliquity cycle (red line). I. The precession cycle extracted from the Induration record (black line) and the Hilbert transform of the precession cycle (red line).
3.7 The phase relationships
The 405-kyr eccentricity cycle was extracted from the induration record and extracted from the Hilbert transform of the 100-kyr eccentricity cycle of the induration record. The 100-kyr eccentricity cycle was extracted from the induration record and extracted from the Hilbert transform of the precession cycle of the induration record. The relationship between the 405-kyr eccentricity cycle directly extracted from the record and the one obtained from the Hilbert transform of the 100-kyr cycle shows an anti-phased relationship. The relationship between the 100-kyr eccentricity cycle directly extracted from the induration record and that obtained from the Hilbert transform of the precession cycle demonstrates a generally anti-phased relationship as well (see Figure 8).
Code
ylims <-c(min(kosov_time[, 2])-5,max(kosov_time[,2])+5)ylims <-rev(ylims/1000)layout.matrix <-matrix(c(1,2),nrow =1,ncol =2 ,byrow =TRUE)graphics::layout(mat = layout.matrix,heights =c(1,1),# Heights of the two rowswidths =c(1,1))par(mar =c(4, 4, 4, 2))plot(x = induration_filt_time_405[, 2],y = induration_filt_time_405[, 1]/1000,type ="l",ylim =ylims,xlab ="405-kyr ecc",col="black",ylab ="Age (Ma)",main="A")lines(x = (induration_filt_time_110_hilb_405[, 2]-mean(induration_filt_time_110_hilb_405[, 2]))*2+mean(induration_filt_time_405[,2]), y = induration_filt_time_110_hilb_405[, 1]/1000, col ="red")plot(x = induration_filt_time_110[, 2],y = induration_filt_time_110[, 1]/1000,type ="l",ylim =ylims,xlab ="100-kyr ecc",col="black",ylab ="Age (Ma)",main="B")lines(x = (induration_filt_time_prec_hilb_110[, 2]-mean(induration_filt_time_prec_hilb_110[, 2]))*2+mean(induration_filt_time_110[,2]), y = induration_filt_time_prec_hilb_110[, 1]/1000, col ="red")
Figure 8. Astronomical cycles extracted from the tuned induration record. A. The black line is a 405-kyr eccentricity cycle extracted from the Induration record. The red line is the 405-kyr eccentricity cycle extracted from the Hilbert transform of the 100-kyr eccentricity cycle of the Induration record. B. The black line is a 100-kyr eccentricity cycle extracted from the Induration record. The red line is the 100-kyr eccentricity cycle extracted from the Hilbert transform of the precession cycle of the Induration record.
3.8 Rate of change in the δ13Ccarb curve
The rate of change (‰/kyr) of the δ13Ccarb curve was calculated from the LOWESS regression applied to the δ¹³Ccarb record. The regression was applied to minimise the influence of short-term local noise (Figure 9). The rate of change (‰/kyr) reaches a maximum of 0.09 (‰/kyr) during the onset of the Lau Biogeochemical Event. The rate of change (‰/kyr) record contains a strong imprint of the 100-kyr eccentricity cycle, which was subsequently extracted from the record and plotted next to the 100-kyr eccentricity cycle extracted from the induration record. During the Lau, the rate of change (‰/kyr) recorded its largest amplitude variations. The 100-kyr eccentricity cycle extracted from the rate of change (‰/kyr) record is generally anti-phased with the 100-kyr eccentricity cycle extracted from the induration record.
Figure 9. Rate of change (‰/kyr) of the δ13Ccarb record. A Stages. B. Event zone subdivision. C. The Induration record. D. The δ13Ccarb record (light green), and the LOWESS smoothed curve (dark green) E. The δ13Ccarb rate of change record (‰/kyr) and the 100-kyr eccentricity cycle extracted from said record F. The 100-kyr eccentricity cycle extracted from the δ13Ccarb rate of change record (‰/kyr) (dark green) and the 100-kyr eccentricity cycle extracted from the induration record.
3.9 The lag-1 sea-level curve
The lag-1 function was used to calculate the lag-1 autocorrelation coefficient using a windowed analysis Monte-Carlo analysis. The resulting curve serves as a proxy for sea-level changes (Li et al., 2018). The lag-1 does show a few trends (see Figure 10). Before the Lau excursion, the curve is characterized by relatively high values (high sea-level) superimposed by some low-amplitude fluctuations. At the onset of the Lau event, a minor rise is observed, followed by a large drop towards a minimum (low sea-level) during the middle of the Lau δ13Ccarb excursion. Towards the top of the studied succession, the curve returns to higher and more regular values (high stable sea-level). The lag-1 record contains the imprint of the 405-kyr eccentricity cycle, which was extracted and compared to the 405-kyr eccentricity cycle extracted from the induration record. The 405-kyr eccentricity cycle extracted from the rate lag-1 record is generally anti-phased with the 405-kyr eccentricity cycle extracted from the induration record.
This code takes a long time to run so one can also download the results from github induration_time_lag_1 <- lag_1( data = induration_time_2, n_sim = 1000, run_multicore = TRUE, win_max = 300, win_min = 50, verbose = TRUE ) write.csv(induration_time_lag_1,“induration_time_lag_1.csv”)
Figure 10. Rate of change of the δ13Ccarb record. A Stages. B. Event zone subdivision. C. The Induration record. D. The δ13Ccarb record (light green) and the LOWESS smoothed curve (dark green) E. The lag-1 sea-level curve (red) and the 405-kyr eccentricity cycle extracted from said record (brown). F. The 405-kyr eccentricity cycle extracted from the lag-1 curve (brown) and the 405-kyr eccentricity cycle extracted from the Induration record (blue).
4. Discussion
4.1 The astrochronology and the imprint of astronomical cycles
The establishment of the the astrochronological framework for the Kosov quarry section is a new step into the quantification of the tempo of Silurian environmental change to show complete the Silurian astronomical timescale. The new anchored astrochronology provides the first near complete high-resolution temporal calibration which allows us to examine the imprint of astronomical cycles and their role in modulating the Lau Biogeochemical Event. The anchored astrochronological age model was used to assign ages, durations, and associated uncertainties to chrono- and geochronologic units, lithological units, and events in the Kosov quarry section, as tabulated in Table 1. The Lau Biogeochemical Event started at 424.06 ± 0.55 Ma and lasted 0.87 ± 0.03 Myr; in contrast, the LKB Event was as short as 30 ± 10 kyr. The durations and associated uncertainties for the Lau Event and its subdivisions are the first astrochronological constraints available for a near-complete record spanning the Lau Event. The durations and associated uncertainties are therefore the best constrained reference values to date. A previous estimate for the onset of the Lau Event is based on the dating of a hiatus in the Cellon section, spanning 423.83 ± 0.55 Ma to 424.19 ± 0.55 Ma, with the Lau Biogeochemical Event starting within the hiatus (Arts et al., 2024). Our age for the onset of the Lau Biogeochemical Event (424.06 ± 0.55 Ma) falls within this range, supporting alignment with previous results. The uncertainties assigned to durations in the record are around ~15%. Since only one proxy was used, a multi-proxy study might even further reduce the uncertainty in the underlying astrochronology of the Kosov quarry section.
The CWT in the time domain contains spectral peaks which can be linked to the 8.5-kyr half precession, 16-kyr precession, 28-kyr obliquity,50-kyr eccentricity, 100-kyr eccentricity, 200-kyr eccentricity and the 405-kyr eccentricity cycle (see Figure 7). The observation of a whole suite of spectral peaks corresponding to known astronomical cycles in the CWT scalogram validates the underlying age model. The periods of the obliquity and precession cycles are at the shorter end of the typically expected durations for the Silurian, yet remain within the current uncertainty range for the estimated durations of these cycles during the Silurian (Waltham, 2015; Farhat et al., 2022; Wu et al., 2024). The spectral peak of ~55 kyr was observed in the average CWT scalogram. This spectral peak can be either the 55 or 54-kyr eccentricity cycle or the 52.8-kyr obliquity cycle (Laskar et al., 2004; Laskar, 2020, p. 20) (see Figure 7). The spectral power of the 55-kyr cycle is highly variable, indicating that the 55-kyr cycle might be of a transient nature unrelated to astronomical forcing. A 200-kyr cycle can also be identified in the record. Given the strong imprint of the 100-kyr and 405-kyr eccentricity cycles, it is reasonable to interpret the 200-kyr cycle as the seldom observed 200-kyr eccentricity cycle (Hilgen et al. (2020)).
Comparable periodicities have been identified in other Silurian successions, such as the Gotland and Yangtze Platform sections, where 100- and 405-kyr cycles dominate δ¹³C and lithological rhythms [citations]. The slight shortening of precession and obliquity cycles in the Kosov record likely reflects uncertainties in tidal dissipation rates and rotational parameters during the mid-Silurian, consistent with model predictions by Waltham (2015) and Wu et al. (2024).
Since amplitude modulating cycles apply unidirectional to the amplitude of the signal, irrespective of whether the proxy exhibits a positive or negative response to insolation forcing it is possible to infer the phase relationship between astronomical cycles extracted from the proxy record and its relationship with the true phase of the astronomical cycle, as well as its relationship with changes in insolation (Hinnov, 2000; Laskar et al., 2004). In the induration record, the precession and short eccentricity amplitude modulations display an anti-phase relationship (see Figure 8). This indicates that during isolation minima, induration values are at a maximum. The anti-phase behavior between induration and insolation suggests that induration intensity in the Kosov quarry reflects carbonate productivity and terrigenous dilution. During insolation minima, reduced runoff and enhanced carbonate precipitation led to higher induration values, while during insolation maxima, seasonality increases, leading to higher runoff, increased detrital products, and decreased carbonate productivity resulting in low induration values (Mutterlose and Ruffell, 1999; Martinez, 2018).
The consistent imprint of eccentricity and precession cycles in both the induration and δ¹³Ccarb records underscores the pervasive sole that astronomical forcing played in pacing of depositonal setting and carbon cycle during the Lau Event. The established relationship between isolation and induration can also be used to understand the imprint of astronomical cycles in the rate of change ‰/kyr) of the δ¹³Ccarb record and the lag-1 sea-level curve. Since the 100-kyr eccentricity cycle extracted from the rate of change (‰/kyr) record is generally anti-phased with the 100-kyr eccentricity cycle extracted from the induration record, this indicates that the rate of change (‰/kyr) is the highest during eccentricity maxima. The 405-kyr eccentricity cycle extracted from the rate lag-1 record is generally anti-phased with the 405-kyr eccentricity cycle extracted from the induration record, indicating that sea-level is the highest during eccentricity maxima.
4.2 An astronomically paced carbon cycle during the Lau Biogeochemical Event
Quantifying the rate of δ¹³Ccarb change (‰/kyr) provides direct insight into the sensitivity and feedback timescales of the Silurian carbon cycle. Until now, the tempo of major Paleozoic carbon isotope excursions remained largely unconstrained, hindering efforts to compare their dynamics with younger Mesozoic and Cenozoic analogues. The new astrochronologically anchored age model allows for the first robust estimation of δ¹³Ccarb change rates (‰/kyr) across the Lau Event, enabling temporal comparison with other Phanerozoic perturbations (Cramer and Jarvis, 2020).
The new age model shows that the rate of change in the δ¹³Ccarb record, reached maximum values of 0.09‰ per kyr for the Lau Event. These rates were derived from a LOWESS regression of the δ¹³Ccarb curve to minimize the influence of short-term local noise. As such, the values represent a robust regional and potentially global benchmark for the Lau event. The high rate of change in the δ¹³Ccarb record (0.09‰ per kyr) during the Lau Event suggests that the Silurian carbon cycle was highly dynamic during the Lau Event especially so compared to the δ¹³Ccarb rate of change curve calculated from the reference curve (Cramer and Jarvis, 2020) (see Figure 11). Considering the underlying age-model and smoothing our estimates for the δ¹³Ccarb rate of change (‰/kyr) presents a conservative rate for the Lau Event. Since it is now possible to constrain the rate of change, it will be possible to model the Lau Biogeochemical event in the domain. This will enable us to examine whether steady state depth domain models, such as those of Frýda et al. (2020). Incorporating the observed δ¹³C change rates into box or Earth system models will help us constrain the response time of the Silurian carbon reservoir and test whether steady-state carbon cycle models remain valid.
The synchronicity between high δ¹³Ccarb change rates and eccentricity maxima suggests a strong astronomically paced modulation of carbon cycle feedback mechanisms. The 100-kyr eccentricity imprint in the δ¹³Ccarb rate-of-change record during the Lau Event indicates that the carbon cycle was at least partly paced by this cycle. The phase relationship shows that the highest rates of change occurred during eccentricity maxima. This suggests that eccentricity maxima resulted in enhanced organic carbon burial during the Lau Event. The observed phase relationship aligns with a scenario in which eccentricity maxima increased the amplitude of climate precession, resulting in enhanced seasonality (Rossignol-Strick, 1985; Van Os et al., 1994; Sachs and Repeta, 1999; Wang, 2009; Gambacorta et al., 2018). The resulting increase in seasonality intensified the hydrological cycle, increasing physical and chemical weathering, supplying fine-grained sediments and nutrients through enhanced (monsoonal) rainfall. The increased fluvial discharge promoted seasonal productivity blooms, facilitating the accumulation of organic-rich deposits. In contrast, during eccentricity minima, the reduced precession amplitude weakened monsoon intensity, lowered nutrient delivery, and re-established more oxic conditions lower carbon burial and thus lowering the δ¹³Ccarb rate-of-change. The enhanced organic carbon burial during the eccentricity maxima would have led to a drawdown of atmospheric CO₂, promoting global cooling and potentially contributing to glacial expansion during the Lau Biogeochemical Evant (Frýda et al., 2021b). The observed coupled coupling between orbital forcing, carbon burial, and climatic feedbacks thus can thus provide a causal linkg between eccentricity forcing, oceanic anoxia and the onset of the Mid-Ludfordian Glaciation Frýda et al. (2021b)].
Figure 11. The GTS 2020 its δ13Ccarb curve and the δ13Ccarb rate of change curve calculated from the δ13Ccarb curve, which is compared to the δ13Ccarb record and δ13Ccarb rate of change record from the Kosov Quarry section.A.Silurian timescale GTS 2020 (shifted up by 0.25 Myrs). B. Silurian δ13Ccarb curve from the GTS 2020 (shifted up by 0.25 Myrs) overlain by the δ13Ccarb curve of the Kosov Quarry section. C . Silurian δ13Ccarb rate of change curve calculated based on the δ13Ccarb curve from the GTS 2020 (shifted up by 0.25 Myrs) overlain by the δ13Ccarb rate of change curve based on the one of the Kosov Quarry section. D. The δ13Ccarb rate of change record (‰/kyr) and the 100-kyr eccentricity cycle extracted from said record.
4.3 Astronomical pacing the lag-1 sea-level curve
The lag-1 curve used as a sea-level proxy aligns with previous sea-level interpretations made for the record of the Kosov quarry section (Manda and Kříž, 2006; Frýda et al., 2020, 2021b) (see figure 10). The sea-level regression inferred from the Kosov record is consistent with global evidence for glacio-eustatic sea-level drawdown during the Lau Event, including sedimentary and isotopic signatures from Laurentia, Baltica, and Gondwana (Manda and Kříž, 2006; Jeppsson et al., 2007; Eriksson and Calner, 2008; Kiipli et al., 2010; Trotter et al., 2016; Frýda et al., 2021b). Our new results when placed ina global context underscored the gloabal synchronicity climatic perturbation assocaited with the Lau Biogeochemical Event and brings forth additaionl supports that a short-lived icehouse episode occured within an otherwise greenhouse like world of the late Silurian.
The 405-kyr eccentricity cycle extracted from the rate lag-1 record is generally anti-phased with the 405-kyr eccentricity cycle extracted from the induration record, indicating that sea level is highest during eccentricity maxima (Figure 10). The relationship between astronomical forcing and sea-level is much akin to that observed when sea-level was under a glacio-eustatic regime (Lourens and Hilgen, 1997; Li et al., 2018). The transition from relatively stable low-amplitude variability into pronounced eccentricity-paced oscillations during the Lau Event highlights the coupling between astronomical forcing, glacio-eustatic change, and the global carbon cycle during the Lau Event. The increased amplitude of sea-level variability indicates that astronomically paced glacio-eustatic fluctuations were amplified by the climatic changes associated with the Lau Biogeochemical Event that culmulated in the Mid-Ludfordian Glaciation (Frýda et al., 2021b).
5. Conclusions
The new astrochronology from the Kosov Quarry provides the first high-resolution temporal framework for the Lau Biogeochemical Event. Cyclostratigraphic analysis demonstrates that the Lau Biogeochemical Event started at 424.03 ± 0.55 Ma and lasted 0.87 ± 0.09 Myr, while the associated LKB Event was much shorter at 30 ± 10 kyr. Astronomical forcing is strongly imprinted across multiple proxies. The induration record is imprinted by astronomical cycles ranging from the 9-kyr half-precession to the 405-kyr eccentricity cycle. The δ¹³Ccarb rate-of-change record contains a strong imprint of the 100-kyr eccentricity cycle, whereas the lag-1 sea-level curve displays a moderate imprint by the 405-kyr eccentricity cycle. The δ¹³Ccarb rate-of-change curve further documents maximum values of 0.09 ‰/kyr during the onset of the Lau, highlighting the highly dynamic character of the Silurian carbon cycle. The imprint of the 100-kyr eccentricity cycle in the δ¹³Ccarb rate-of-change also indicates that this cycle acted as a pacemaker of carbon cycle variability and carbon burial during the Lau Event. The strong imprint of the 405-kyr eccentricity cycle in the lag-1 sea-level and the inferred phase relationship indicate that sea-level was subjected to an astronomically paced glacio-eustatic regime. Collectively, these findings demonstrate that the Lau Biogeochemical Event represents a climatically amplified eccentricity-paced perturbation of the Silurian carbon cycle. The astrochronologically constrained rate of δ¹³C change, coupled with coeval glacio-eustatic oscillations, highlights the sensitivity of the mid-Paleozoic Earth system to orbital forcing. These results suggest that even during intervals of high baseline warmth, astronomically paced feedbacks between climate dynamics, organic carbon burial, and ice-sheet growth could operate effectively, shaping both carbon cycling and sea-level evolution on eccentricity timescales.
Figure SI 1. Induration proxy record in the time domain, including the wavelet scalogram of the induration record and astronomical cycles extracted from said record. A Stages. B. Event zone subdivision C. The δ13Ccarb record. D. The Induration record. E. Wavelet scalogram of the induration record with the average spectral power on top. The black vertical lines in the wavelet scalograms are durations of known astronomical cycles. From left to right, these cycles are the 8.5-kyr half precession, ~16-kyr precession, 28-kyr obliquity,50-kyr eccentricity, 100-kyr eccentricity, 200-kyr eccentricity and the 405-kyr eccentricity cycle. F. The black line is a 405-kyr eccentricity cycle extracted from the Induration record. The red line is the 405-kyr eccentricity cycle extracted from the Hilbert transform of the 100-kyr eccentricity cycle of the Induration record. G. The 100-kyr eccentricity cycle was extracted from the Induration record (black line) and the Hilbert transform of the 100-kyr eccentricity cycle (red line). H.) Obliquity cycle extracted from the Induration record (black line), and the Hilbert transform of the obliquity cycle (red line). I. The precession cycle extracted from the Induration record (black line) and the Hilbert transform of the precession cycle (red line).
6. References
Arts, M., Corradini, C., Pondrelli, M., Pas, D., and Da Silva, A.C., 2024, Age and orbital forcing in the upper Silurian Cellon section (Carnic Alps, Austria) uncovered using the WaverideR R package: Frontiers in Earth Science, v. 12, p. 1–24, doi:10.3389/feart.2024.1357751.
Cramer, B.D., Schmitz, M.D., Huff, W.D., and Bergström, S.M., 2015, High-precision u–pb zircon age constraints on the duration of rapid biogeochemical events during the ludlow epoch (Silurian period): Journal of the Geological Society, v. 172, p. 157–160, doi:10.1144/jgs2014-094.
Eriksson, M.J., and Calner, M., 2008, A sequence stratigraphical model for the Late Ludfordian (Silurian) of Gotland, Sweden: Implications for timing between changes in sea level, palaeoecology, and the global carbon cycle: Facies, v. 54, p. 253–276, doi:10.1007/s10347-007-0128-y.
Farhat, M., Auclair-Desrotour, P., Boué, G., and Laskar, J., 2022, The resonant tidal evolution of the Earth-Moon distance: Astronomy & Astrophysics, v. 665, p. L1.
Fatka, O., and Mergl, M., 2009, The “microcontinent”Perunica: Status and story 15 years after conception: Geological Society, London, Special Publications, v. 325, p. 65–101, doi:10.1144/SP325.4.
Frýda, J., Lehnert, O., Frýdová, B., Farkaš, J., and Kubajko, M., 2020, Carbon and sulfur cycling during the mid-Ludfordian anomaly and the linkage with the late Silurian Lau/Kozlowskii Bioevent: Palaeogeography, Palaeoclimatology, Palaeoecology, v. 564, p. 110152, doi:10.1016/j.palaeo.2020.110152.
Frýda, J., Lehnert, O., Joachimski, M.M., Männik, P., Kubajko, M., Mergl, M., Farkaš, J., and Frýdová, B., 2021b, The Mid-Ludfordian (late Silurian) Glaciation: A link with global changes in ocean chemistry and ecosystem overturns: Earth-Science Reviews, v. 220, p. 103652, doi:10.1016/j.earscirev.2021.103652.
Frýda, J., Lehnert, O., Joachimski, M.M., Männik, P., Kubajko, M., Mergl, M., Farkaš, J., and Frýdová, B., 2021a, The Mid-Ludfordian (late Silurian) Glaciation: A link with global changes in ocean chemistry and ecosystem overturns: Earth-Science Reviews, v. 220, p. 103652, doi:10.1016/j.earscirev.2021.103652.
Frýda, J., and Manda, Š., 2013, A long-lasting steady period of isotopically heavy carbon in the late Silurian ocean: Evolution of the δ13 record and its significance for an integrated δ13, graptolite and conodont stratigraphy: Bulletin of Geosciences, p. 463–482, doi:10.3140/bull.geosci.1436.
Gambacorta, G., Menichetti, E., Trincianti, E., and Torricelli, S., 2018, Orbital control on cyclical primary productivity and benthic anoxia: Astronomical tuning of the Telychian Stage (Early Silurian): Palaeogeography, Palaeoclimatology, Palaeoecology, v. 495, p. 152–162, doi:10.1016/j.palaeo.2018.01.003.
Hilgen, F., Zeeden, C., and Laskar, J., 2020, Paleoclimate records reveal elusive 200-kyr eccentricity cycle for the first time: Global and Planetary Change, v. 194, p. 103296, doi:10.1016/j.gloplacha.2020.103296.
Hinnov, L.A., 2000, New Perspectives on Orbitally Forced Stratigraphy: Annual Review of Earth and Planetary Sciences, v. 28, p. 419–475, doi:10.1177/0741713604268894.
HORNÝ, R., 1955, Base vrstev kopaninských eß 1 na jihozápadním okraji vulkanické facie (Kosov u Berouna): Věstník Ústředního ústavu geologického, v. 30, p. 81–86.
Jeppsson, L., Talent, J.A., Mawson, R., Simpson, A.J., Andrew, A.S., Calner, M., Whitford, D.J., Trotter, J.A., Sandström, O., and Caldon, H.J., 2007, High-resolution Late Silurian correlations between Gotland, Sweden, and the Broken River region, NE Australia: Lithologies, conodonts and isotopes: Palaeogeography, Palaeoclimatology, Palaeoecology, v. 245, p. 115–137, doi:10.1016/j.palaeo.2006.02.032.
Kiipli, T., Kiipli, E., and Kaljo, D., 2010, Silurian sea level variations estimated using SiO2/Al2O 3 and K2O/A12O3 ratios in the Priekule drill core section, latvia: Bollettino della Societa Paleontologica Italiana, v. 49, p. 55–63.
Kříž, J., 1992, Silurian field excursions: Prague Basin (Barrandian), Bohemia: National Museum of Wales, v. 13.
Kříž, J., 1991, The Silurian of the Prague Basin (Bohemia)–tectonic, eustatic and volcanic controls on facies and faunal development: Special Papers in Palaeontology, v. 44, p. 179–203.
Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, a.C.M., and Levrard, B., 2004, A long-term numerical solution for the insolation quantities of the Earth: Astronomy and Astrophysics, v. 428, p. 261–285, doi:10.1051/0004-6361:20041335.
Lehnert, O., Frýda, J., Buggisch, W., Munnecke, A., Nützel, A., Křiž, J., and Manda, S., 2007, δ13C records across the late Silurian Lau event: New data from middle palaeo-latitudes of northern peri-Gondwana (Prague Basin, Czech Republic): Palaeogeography, Palaeoclimatology, Palaeoecology, v. 245, p. 227–244, doi:10.1016/j.palaeo.2006.02.022.
Li, M., Hinnov, L.A., Huang, C., and Ogg, J.G., 2018, Sedimentary noise and sea levels linked to land-ocean water exchange and obliquity forcing: Nature Communications, v. 9, doi:10.1038/s41467-018-03454-y.
Lourens, L.J., and Hilgen, F.J., 1997, Long-periodic variations in the earth’s obliquity and their relation to third-order eustatic cycles and late Neogene glaciations: Quaternary International, v. 40, p. 43–52, doi:10.1016/S1040-6182(96)00060-2.
Manda, Š., and Kříž, J., 2006, Environmental and biotic changes in subtropical isolated carbonate platforms during the Late Silurian Kozlowskii Event, Prague Basin: GFF, v. 128, p. 161–168.
Martinez, M., 2018, Mechanisms of Preservation of the Eccentricity and Longer-term Milankovitch Cycles in Detrital Supply and Carbonate Production in Hemipelagic Marl-Limestone Alternations: Elsevier Ltd, v. 3, 189–218 p., doi:10.1016/bs.sats.2018.08.002.
Melchin, M.J., Sadler, P.M., and Cramer, B.D., 2020, The Silurian Period (F. M. Gradstein, J. G. Ogg, M. D. Schmitz, & G. M. Ogg, Eds.): BV, 695–732 p., doi:10.1016/B978-0-12-824360-2.00021-8.
Meyers, S.R., 2019, Cyclostratigraphy and the problem of astrochronologic testing: Earth-Science Reviews, v. 190, p. 190–223, doi:10.1016/j.earscirev.2018.11.015.
Meyers, S.R., 2015, The evaluation of eccentricity-related amplitude modulation and bundling in paleoclimate data: An inverse approach for astrochronologic testing and time scale optimization: Paleoceanography, v. 30, p. 1625–1640, doi:10.1002/2015PA002850.
Mutterlose, J., and Ruffell, A., 1999, Milankovitch-scale palaeoclimate changes in pale-dark bedding rhythms from the Early Cretaceous (Hauterivian and Barremian) of eastern England and northern Germany: Palaeogeography, Palaeoclimatology, Palaeoecology, v. 154, p. 133–160, doi:10.1016/S0031-0182(99)00107-8.
Sachs, J.P., and Repeta, D.J., 1999, Oligotrophy and nitrogen fixation during eastern Mediterranean sapropel events: Science, v. 286, p. 2485–2488.
Scotese, C.R., 2021, An atlas of phanerozoic paleogeographic maps: The seas come in and the seas go out: Annual Review of Earth and Planetary Sciences, v. 49, p. 679–728, doi:10.1146/annurev-earth-081320-064052.
Štorch, P., 1995, Biotic crises and post-crisis recoveries recorded by Silurian planktonic graptolite faunas of the Barrandian area (Czech Republic): Geolines, v. 3, p. 59–70.
Tasáryová, Z., Schnabl, P., Čížková, K., Pruner, P., Janoušek, V., Rapprich, V., Štorch, P., Manda, Š., Frýda, J., and Trubač, J., 2014, Gorstian palaeoposition and geotectonic setting of Suchomasty Volcanic Centre (Silurian, Prague Basin, Teplá-Barrandian Unit, Bohemian Massif): Gff, v. 136, p. 262–265, doi:10.1080/11035897.2013.879735.
Trotter, J.A., Williams, I.S., Barnes, C.R., Männik, P., and Simpson, A.J., 2016, New conodont δ18O records of Silurian climate change: Implications for environmental and biological events: Palaeogeography, Palaeoclimatology, Palaeoecology, v. 443, p. 34–48, doi:10.1016/j.palaeo.2015.11.011.
Van Os, B., Lourens, L., Hilgen, F., De Lange, G., and Beaufort, L., 1994, The formation of Pliocene sapropels and carbonate cycles in the Mediterranean: Diagenesis, dilution, and productivity: Paleoceanography, v. 9, p. 601–617.
Waltham, D., 2015, Milankovitch Period Uncertainties and Their Impact On Cyclostratigraphy: Journal of Sedimentary Research, v. 85, p. 990–998, doi:10.2110/jsr.2015.66.
Wang, P., 2009, Global monsoon in a geological perspective: Chinese Science Bulletin, v. 54, p. 1113–1136.
Wu, Y., Malinverno, A., Meyers, S.R., and Hinnov, L.A., 2024, A 650-Myr history of Earth’s axial precession frequency and the evolution of the Earth-Moon system derived from cyclostratigraphy: Science Advances, v. 10, p. eado2412.
